function r = hydro(Ti, TIME_STEPS)
  % First order hydro
  %  eliminating the u-dot-dot terms 
  %  fsolve search method
  close all;
  fig = figure; 
  load('cm', 'mycmap');
  set(fig, 'Colormap', mycmap);

  HYDROLIZE = false;

  GRID_SIZE = 8;
  DX = 0.01;
  Dt = 0.001;
  VECTOR_SUBSAMPLING = 4;
  [gridx, gridy] = meshgrid([1:VECTOR_SUBSAMPLING:GRID_SIZE], [1:VECTOR_SUBSAMPLING:GRID_SIZE]);


  VIDEO = false;
  FPS = 30;
  SECONDS_PER_LIGHTYEAR = 15;
  if (VIDEO)
    vt = VideoWriter('v_t.avi', 'Archival');
    open(vt);
  end

  ZERO_HYDRO = false;



  A = sqrt(3)*pi/18 - log(3)/2;
  c = -1;
  h = (A + 1);
  l = (2 * A);
  
  ux = zeros(GRID_SIZE, GRID_SIZE);
  uy = zeros(GRID_SIZE, GRID_SIZE);
  ut = ones(GRID_SIZE, GRID_SIZE);
  T  = zeros(GRID_SIZE, GRID_SIZE);
  PI = zeros(GRID_SIZE, GRID_SIZE);
  B = zeros(GRID_SIZE, GRID_SIZE);
  
  ux0 = zeros(GRID_SIZE, GRID_SIZE);
  uy0 = zeros(GRID_SIZE, GRID_SIZE);
  T0  = zeros(GRID_SIZE, GRID_SIZE);
  PI0 = zeros(GRID_SIZE, GRID_SIZE);
  B0 = zeros(GRID_SIZE, GRID_SIZE);
  
  DT  = zeros(GRID_SIZE, GRID_SIZE);
  Dux = zeros(GRID_SIZE, GRID_SIZE);
  Duy = zeros(GRID_SIZE, GRID_SIZE);
  DPI = zeros(GRID_SIZE, GRID_SIZE);
  DB = zeros(GRID_SIZE, GRID_SIZE);
  
  dtT  = zeros(GRID_SIZE, GRID_SIZE); 
  dtux = zeros(GRID_SIZE, GRID_SIZE); 
  dtuy = zeros(GRID_SIZE, GRID_SIZE); 

  dxT = zeros(GRID_SIZE, GRID_SIZE); 
  dyT = zeros(GRID_SIZE, GRID_SIZE); 
  dxux = zeros(GRID_SIZE, GRID_SIZE); 
  dyux = zeros(GRID_SIZE, GRID_SIZE); 
  dxuy = zeros(GRID_SIZE, GRID_SIZE); 
  dyuy = zeros(GRID_SIZE, GRID_SIZE);
  dxut = zeros(GRID_SIZE, GRID_SIZE); 
  dyut = zeros(GRID_SIZE, GRID_SIZE);
  dxPI = zeros(GRID_SIZE, GRID_SIZE);
  dyPI = zeros(GRID_SIZE, GRID_SIZE);
  dxB = zeros(GRID_SIZE, GRID_SIZE);
  dyB = zeros(GRID_SIZE, GRID_SIZE);

  
  
  
  dxdxT = zeros(GRID_SIZE, GRID_SIZE); 
  dxdyT = zeros(GRID_SIZE, GRID_SIZE); 
  dydyT = zeros(GRID_SIZE, GRID_SIZE); 

  dxdxux = zeros(GRID_SIZE, GRID_SIZE); 
  dxdyux = zeros(GRID_SIZE, GRID_SIZE); 
  dydyux = zeros(GRID_SIZE, GRID_SIZE); 

  dxdxuy = zeros(GRID_SIZE, GRID_SIZE); 
  dxdyuy = zeros(GRID_SIZE, GRID_SIZE); 
  dydyuy = zeros(GRID_SIZE, GRID_SIZE); 



  sumT = 0;
  % initial conditions for vx vy come here
  % vx
  % vy
  % T
  TEMP_OFFSET =0.7;
  cx = 0.5 * GRID_SIZE + 0.5;
  cy = 0.5 * GRID_SIZE + 0.5;
  for x = 1:GRID_SIZE
    for y = 1:GRID_SIZE
      % T(x, y) = 0.00001 * exp(-((x-cx)^2 + (y-cy)^2)/12);
      % T(x,y) = 0.2*sin(2*pi*x/GRID_SIZE);
	  % T(x, y) = 0.001
	   T(x, y) = 0.1 * sin(2*pi*(x/GRID_SIZE));
    end
  end
  % T(:,:) = Ti(:,:);
  T = T + TEMP_OFFSET;
   %for x = 1:GRID_SIZE
    %  for y = 1:GRID_SIZE
       %T(x,y) = T(x,y) +  0.01 * rand();
    %   ux(x,y) = 0.00001 * sin(2*pi*(x/GRID_SIZE));
       %vy(x,y) = 0.001 * rand();
    %  end
   %end

  %vx = 0.000000 * ones(GRID_SIZE, GRID_SIZE);;
  %vy = 0.2 * ones(GRID_SIZE, GRID_SIZE);;
  
  % ****************************************************************
  %  hydrolize the initial conditions
  % ****************************************************************
  if (HYDROLIZE)
    level = 4;
    k = zeros(GRID_SIZE+1, GRID_SIZE+1);  
    [cx, cy] = meshgrid(0:GRID_SIZE-1,0:GRID_SIZE-1);
    [kx, ky] = meshgrid(-GRID_SIZE/2:GRID_SIZE/2);
    for kx = -GRID_SIZE/level:GRID_SIZE/level
      for ky = -GRID_SIZE/level:GRID_SIZE/level
        ft = T .* exp(-2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
        k(kx+GRID_SIZE/2+1, ky+GRID_SIZE/2+1) = sum(sum(ft));
      end
    end
    T = zeros(GRID_SIZE:GRID_SIZE);
    for kx = -GRID_SIZE/level:GRID_SIZE/level
      for ky = -GRID_SIZE/level:GRID_SIZE/level
        T = T + k(kx+GRID_SIZE/2+1,ky+GRID_SIZE/2+1) * exp(2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
        % vx = vx + 20*(rand()-0.5)*exp(2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
        % vy = vy + 20*(rand()-0.5)*exp(2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);

      end
    end
    T = real(T / GRID_SIZE / GRID_SIZE);
    %ux = real(ux / GRID_SIZE / GRID_SIZE);
    %uy = real(uy / GRID_SIZE / GRID_SIZE);
  end



  % ****************************************************************
  %  compute derivative kernels
  % ****************************************************************
  % can upgrade on this -- trefethen pg 6

  d1 = zeros(GRID_SIZE, GRID_SIZE);
  d2 = zeros(GRID_SIZE, GRID_SIZE);
  for x = 1:GRID_SIZE
    for y=1:GRID_SIZE
      s1 = 0;
      s2 = 0;
      for k = -GRID_SIZE/2 : GRID_SIZE/2
        s1 = s1 + 2. * pi * i / GRID_SIZE / GRID_SIZE * k * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
        s2 = s2 + (2. * pi * i / GRID_SIZE * k)^2 / GRID_SIZE * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
      end
      d1(x,y) = real(s1)/DX;
      d2(x,y) = real(s2)/DX/DX;
    end
  end
  h = 2*pi/GRID_SIZE;
  col1 = [0 .5*(-1).^(1:GRID_SIZE-1).*cot((1:GRID_SIZE-1)*h/2)]  * h/DX;
  max(max(d1 - toeplitz(col1, col1([1 GRID_SIZE:-1:2]))))
  col2 = [(-pi^2/(3*h*h)-1/6) .5*(-1).^(0:GRID_SIZE-2).*(csc((1:GRID_SIZE-1)*h/2).^2)]  * h*h/(DX*DX);
  % d2 = toeplitz(col2);
  

  %dxT = dxdx(T);
  %dxT2 = dxdx2(T);
  % image(imagesc(dxT));
  %figure;
  %image(imagesc(dxT2));
  %image((dxT) * 128  + 128); 
  % image((dxT) * 128 + 128); 
  % truesize(fig);
  %fig=figure;
  %truesize(fig);
  %keyboard();
 % image((T-TEMP_OFFSET)*256+128);

  %return

  dxux = dx(ux);
  dyux = dy(ux);
  dxuy = dx(uy);
  dyuy = dy(uy);
  dxT = dx(T);
  dyT = dy(T);
  dxPI = dx(PI);
  dyPI = dy(PI);
  dxB = dx(B);
  dyB = dy(B);
  
  dxdxux = dxdx(ux);
  dxdyux = dx(dy(ux));
  dydyux = dydy(ux);
  dxdxuy = dxdx(uy);
  dxdyuy = dx(dy(uy));
  dydyuy = dydy(uy);
  dxdxT = dxdx(T);
  dxdyT = dx(dy(T));
  dydyT = dydy(T);


  tic;

	
    % solving non-linear equation of type:
    % L(u, T, u-dot, T-dot) + R(u, T) == 0  -- 3 components
    % each of L, R has 3 components 

    % precompute the R functions
    % R1
    a = -1;
    R1 = (1./(2.*ut.^2)).*(B.*dxuy.*ut + B.*dyux.*ut + dxuy.*PI.*ut + dyux.*PI.*ut+3.*dxux.*T.^3.*ut.^3 + 3.*dyuy.*T.^3.*ut.^3 - B.*dyut.*ux - dyut.*PI.*ux + dyB.*ut.*ux + dyPI.*ut.*ux + 3.*dxT.*T.^2.*ut.*ux - 3.*dxut.*T.^3.*ut.^2.*ux - B.*dxuy.*ut.*ux.^2 + 3.*B.*dyux.*ut.*ux.^2 - B.*dyut.*ux.^3 + dyB.*ut.*ux.^3 - (dxut.*(B + PI) - ut.*(dxB + dxPI + 3.*T.^2.*(dyT - dyut.*T.*ut)) - 2.*(B.*dyuy + dxux.*PI).*ut.*ux + (B.*dxut + 2.*dxut.*PI + dxB.*ut).*ux.^2).*uy + (3.*dxuy.*PI.*ut - dyux.*PI.*ut - (2.*B.*dyut + dyut.*PI + dyPI.*ut).*ux).*uy.^2 + ((-dxut).*PI + dxPI.*ut).*uy.^3);
    %a = 1;
    
    R2 = (1./(2.*ut.^3)).*(ut.*((dyB + dyPI + 3.*dxT.*T.^2).*ut.^2 - dxuy.*(B + PI).*ux + (2.*B.*dyux + 2.*dxuy.*PI + 3.*(dxux + dyuy).*T.^3).*ut.^2.*ux + ut.*(-3.*dxut.*T.^3 + dyB.*ut).*ux.^2 - dxuy.*(B + 2.*PI).*ux.^3) - (dyuy.*(B + PI).*ut - 2.*(dxux + dyuy).*PI.*ut.^3 + ((-dxut).*(B + PI) + ut.*(dxB + dxPI + 3.*dyut.*T.^3.*ut - 2.*dxPI.*ut.^2)).*ux + (2.*dxux + dyuy).*(B + 2.*PI).*ut.*ux.^2 + ((-dxut).*(B + 2.*PI) + (dxB + 2.*dxPI).*ut).*ux.^3).*uy + (dyut.*(B + PI) - (dyB + dyPI).*ut + dyPI.*ut.^3 - (2.*B.*dyux + 3.*dxuy.*PI + 4.*dyux.*PI).*ut.*ux + (B.*dyut + 2.*dyut.*PI - (dyB + 2.*dyPI).*ut).*ux.^2).*uy.^2 + (-3.*dyuy.*PI.*ut + dxut.*PI.*ux - dxPI.*ut.*ux).*uy.^3 + (dyut.*PI - dyPI.*ut).*uy.^4);
    %R2 = R2 + a ..* (    dxdx(ux...*dxux + uy...*dyux) + dydy(ux...*dxux + uy...*dyux)  ) + a ..* (    dxdx(ux...*dxux + uy...*dyux) + dx(dy(ux...*dxuy + uy...*dyuy))  );
    %R2 = R2 + a ..* ( ux...*dx(dxux + dyuy) + uy...*dy(dxux + dyuy));

    R3 = (1./(2.*ut.^3)).*(3.*dyT.*T.^2.*ut.^3 - B.*dxux.*ut.*ux - dxux.*PI.*ut.*ux +2.*B.*dxux.*ut.^3.*ux + 2.*B.*dyuy.*ut.^3.*ux + B.*dxut.*ux.^2 + dxut.*PI.*ux.^2 -    3.*B.*dxux.*ut.*ux.^3 + B.*dxut.*ux.^4 +       ((-dyux).*(B + PI).*       ut + (2.*B.*dyux + 2.*dxuy.*PI + 3.*(dxux + dyuy).*T.^3).*       ut.^3 + (dyut.*(B + PI) - (dyB + dyPI).*ut - 3.*dxut.*T.^3.*ut.^2 +          2.*dyB.*ut.^3).*ux -            (4.*B.*dxuy + 3.*B.*dyux + 2.*dxuy.*PI).*ut.*       ux.^2 + (B.*dyut - dyB.*ut).*ux.^3).*    uy - (3.*dyut.*T.^3.*ut.^2 + (2.*B + PI).*       ux.*(dxux.*ut + 2.*dyuy.*ut - dxut.*ux)).*uy.^2 -       (2.*B.*dyux.*ut + dyux.*PI.*ut - 2.*B.*dyut.*ux - dyut.*PI.*ux +       2.*dyB.*ut.*ux + dyPI.*ut.*ux).*uy.^3 +    dxPI.*ut.*(ut - ux).*(ut + ux).*(1 + uy.^2) +       dxB.*ut.*((ut - ux).*(ut + ux).*(1 + ux.^2) - 2.*ux.^2.*uy.^2));
    %R3 = R3 + a ..* (    dx(dy(ux...*dxux + uy...*dyux)) + dydy(ux...*dxuy + uy...*dyuy)  ) + a ..* (    dxdx(ux...*dxuy + uy...*dyuy) + dydy(ux...*dxuy + uy...*dyuy)  );
    %R2 = R2 + a ..* ( ux...*dx(dxux + dyuy) + uy...*dy(dxux + dyuy));

    R4=-((1./(4.*c.*T.*ut.^2)).*(dxuy.*l.*PI.*ut.^2 - dyux.*l.*PI.*ut.^2 +      2.*c.^2.*dxux.*T.^3.*ut.^2 - 2.*c.^2.*dyuy.*T.^3.*ut.^2 +      2.*c.^2.*dxut.*T.^3.*ut.*ux + 2.*dyut.*l.*PI.*ut.^3.*ux -          4.*c.^2.*dxut.*T.^3.*ut.^3.*ux - dyut.*l.*PI.*ut.^5.*ux -      2.*dyux.*l.*PI.*ut.^2.*ux.^2 + 2.*c.^2.*dxux.*T.^3.*ut.^2.*ux.^2 -      2.*c.^2.*dyuy.*T.^3.*ut.^2.*ux.^2 + dyux.*l.*PI.*ut.^4.*ux.^2 +          2.*c.^2.*dxut.*T.^3.*ut.*ux.^3 + dyut.*l.*PI.*ut.^3.*ux.^3 -      dyux.*l.*PI.*ut.^2.*      ux.^4 + (-4.*c.*PI.*T.*ut.^2.*ux +         2.*c.^2.*T.^3.*ut.*(dyut + 2.*dxuy.*ut.*ux + dyut.*ux.^2) -                       dxut.*PI.*ut.*(2.*h.*ux.^2 +            l.*(1 + ux.^2.*(5 + 8.*ux.^2 + 5.*ux.^4 + ux.^6 +                  ut.^4.*(3 + ux.^2) - ut.^2.*(7 + 8.*ux.^2 + 2.*ux.^4)))) +                       PI.*ux.*(dyuy.*ut.^2.*(2.*h + l.*(-2 + ut.^2 - ux.^2)) +            dxux.*(2.*h.*ut.^2 +               l.*(1 + ux.^2.*(5 + 8.*ux.^2 + 5.*ux.^4 + ux.^6 +                     ut.^4.*(3 + ux.^2) - ut.^2.*(7 + 8.*ux.^2 + 2.*ux.^4)))))).*           uy +      PI.*(dyux.*         l.*(-1 + ut.^2 - ux.^2).*(-1 - 4.*ux.^2 + (-4 + ut.^2).*ux.^4 -            ux.^6) +         dxuy.*l.*(ut.^2 + (1 - 5.*ut.^2 + 3.*ut.^4).*            ux.^2 + (6 - 8.*ut.^2 + ut.^4).*ux.^4 +                    (5 - 2.*ut.^2).*ux.^6 + ux.^8) -         dyut.*ut.*ux.*(2.*h +            l.*(1 - 2.*ut.^2 + (6 - 5.*ut.^2 + ut.^4).*ux.^2 + (5 - 2.*ut.^2).*               ux.^4 + ux.^6))).*uy.^2 +          l.*      PI.*(dxut.*ut.*(-1 + ux.^2.*(-1 + ut.^2 - ux.^2).*(3 + 2.*ux.^2)) +         ux.*(dyuy - 2.*dyuy.*ut.^2 +            dyuy.*ux.^2.*(2 - ut.^2 + ux.^2).*(3 - ut.^2 + ux.^2) +                    dxux.*(1 + ux.^2.*(1 - ut.^2 + ux.^2).*(3 + 2.*ux.^2)))).*      uy.^3 + l.*      PI.*(dyux - dyut.*ut.*ux + (dxuy + 3.*dyux - 3.*dxuy.*ut.^2).*ux.^2 +         dyut.*ut.*(-5 + 2.*ut.^2).*ux.^3 -               (dxuy + dyux).*(-5 + 2.*ut.^2).*ux.^4 - 2.*dyut.*ut.*ux.^5 +         2.*(dxuy + dyux).*ux.^6).*uy.^4 +          l.*PI.*      ux.*(dyuy + dyuy.*(5 - 2.*ut.^2).*ux.^2 -         dxut.*ut.*ux.^3 + (dxux + 2.*dyuy).*ux.^4).*uy.^5 +      l.*PI.*ux.^3.*((-dyut).*ut + (dxuy + dyux).*ux).*uy.^6 +      dyuy.*l.*PI.*ux.^3.*uy.^7 +          B.*      l.*(ut.^2.*(1 + ux.^2).*(dxuy + dyut.*ut.*ux.*(2 - ut.^2 + ux.^2) -            dyux.*(1 - (-2 + ut.^2).*ux.^2 + ux.^4)) +               (1 + ux.^2).*(dxux.*ux + dyuy.*ut.^2.*ux.*(-2 + ut.^2 - ux.^2) -            dxut.*ut.*(1 + (3 - 3.*ut.^2 + ut.^4).*ux.^2 + (3 - 2.*ut.^2).*               ux.^4 + ux.^6) +                               dxux.*ux.^3.*(3 + ut.^4 + 3.*ux.^2 + ux.^4 - ut.^2.*(3 + 2.*ux.^2))).*         uy + (dyux +            dyut.*ut.*(-1 + ut.^4).*            ux + (dxuy +               4.*dyux - (dxuy + 2.*dyux).*ut.^2 + (dxuy - dyux).*ut.^4).*                     ux.^2 -            dyut.*ut.*(-2 + ut.^2).^2.*            ux.^3 + (6.*               dyux + (-4 + ut.^2).*(-dxuy + (dxuy + dyux).*ut.^2)).*            ux.^4 + 2.*dyut.*ut.*(-2 + ut.^2).*ux.^5 -                    2.*(dxuy + dyux).*(-2 + ut.^2).*ux.^6 -            dyut.*ut.*ux.^7 + (dxuy + dyux).*ux.^8).*uy.^2 +         ux.*(ux.*((-dxut).*ut + dxux.*ux).*(1 + (3 - 2.*ut.^2).*ux.^2 +               2.*ux.^4) +                               dyuy.*(1 - ut.^4 + (-2 + ut.^2).^2.*ux.^2 - 2.*(-2 + ut.^2).*ux.^4 +               ux.^6)).*uy.^3 +                       ux.*(ux.*(-dxuy + dyux +               2.*dyux.*ut.^2 - (dxuy + dyux).*(-3 + 2.*ut.^2).*ux.^2 +               2.*(dxuy + dyux).*ux.^4) +            dyut.*ut.*(1 - 3.*ux.^2 - 2.*ux.^4 + 2.*ut.^2.*(-1 + ux.^2))).*uy.^4 +         ux.*(ux.*((-dxut).*ut + dxux.*ux).*(-1 + ux.^2) +            dyuy.*(-1 + 3.*ux.^2 + 2.*ux.^4 - 2.*ut.^2.*(-1 + ux.^2))).*uy.^5 +         ux.*((-dyut).*ut + (dxuy + dyux).*ux).*(-1 + ux.^2).*uy.^6 +               dyuy.*ux.*(-1 + ux.^2).*uy.^7)));
	
   

	R5 = -((1./(4.*c.*T.*ut.^2)).*((-B).*dxuy.*l.*ut.^2 + B.*dyux.*l.*ut.^2 -      dxuy.*l.*PI.*ut.^2 + dyux.*l.*PI.*ut.^2 - 2.*c.^2.*dxux.*T.^3.*ut.^2 +      2.*c.^2.*dyuy.*T.^3.*ut.^2 - B.*dyut.*l.*ut.*ux -          dyut.*l.*PI.*ut.*ux + 2.*c.^2.*dxut.*T.^3.*ut.*ux + B.*dxuy.*l.*ux.^2 +      dxuy.*l.*PI.*ux.^2 - B.*dxuy.*l.*ut.^2.*ux.^2 + B.*dyux.*l.*ut.^2.*ux.^2 -      B.*dyut.*l.*ut.*ux.^3 + B.*dxuy.*l.*ux.^4 +          (-4.*B.*c.*T.*ut.^2.*ux +         2.*c.^2.*T.^3.*ut.*(dyut - 2.*dyut.*ut.^2 + 2.*dyux.*ut.*ux) +                       l.*PI.*(dyuy.*             ux + (dxut.*ut - dxux.*ux).*(-ux.^2 + ux.^4 + ux.^6 +               ut.^4.*(-1 + ux.^2) - 2.*ut.^2.*(-1 + ux.^4))) +                       B.*((-dxut).*             ut.*(2.*h.*ux.^2 +               l.*(ut.^4 + ux.^2 + ux.^4 - 2.*ut.^2.*(1 + ux.^2))) +            ux.*(2.*dxux.*h.*ut.^2 + dyuy.*(l + 2.*h.*ut.^2 + l.*ux.^2) +                                       dxux.*l.*(ut.^4 + ux.^2 + ux.^4 - 2.*ut.^2.*(1 + ux.^2))))).*       uy + (ut.^2.*(2.*c.^2.*(-dxux + dyuy).*T.^3 + B.*dxuy.*l.*(-2 + ut.^2) +            l.*PI.*(dyux + dxuy.*(-3 + ut.^2))) -                       ut.*(-2.*c.^2.*dxut.*T.^3 + dyut.*l.*PI.*(4 - 3.*ut.^2 + ut.^4) +            B.*dyut.*(2.*h + l.*(5 - 7.*ut.^2 + 3.*ut.^4))).*ux +                       l.*(5.*B.*dxuy + B.*dyux + 4.*dxuy.*PI +            dyux.*PI - (4.*B.*dxuy + 5.*B.*dyux + 2.*dxuy.*PI + dyux.*PI).*             ut.^2 + (3.*B.*dyux + (-dxuy + dyux).*PI).*ut.^4).*ux.^2 +               dyut.*l.*ut.*(-PI + 3.*B.*(-1 + ut.^2)).*ux.^3 +         l.*(PI.*(dxuy - dyux + 2.*dxuy.*ut.^2) +            B.*(3.*dxuy + dyux - 3.*dyux.*ut.^2)).*ux.^4 + dyut.*l.*PI.*ut.*ux.^5 -               (dxuy + dyux).*l.*PI.*ux.^6).*uy.^2 +          (2.*c.^2.*dyut.*T.^3.*ut +         l.*(dxut.*ut.^3.*(B -               PI.*(-3 + ut.^2)) + (dyuy.*(5.*B +                  4.*PI) - (B.*(dxux + 7.*dyuy) + 3.*(dxux + dyuy).*PI).*                ut.^2 +                         (3.*B.*dyuy + (dxux + dyuy).*PI).*ut.^4).*ux -            dxut.*ut.*(-2 + ut.^2).*(-3.*B - 2.*PI + (B + PI).*ut.^2).*ux.^2 +                    (PI.*(dyuy + dxux.*(-2 + ut.^2).^2) +               B.*(3.*dyuy - 3.*dyuy.*ut.^2 + dxux.*(6 - 5.*ut.^2 + ut.^4))).*             ux.^3 + dxut.*ut.*(-5.*B - 3.*PI + 2.*(B + PI).*ut.^2).*ux.^4 -                    (-5.*B.*dxux - 3.*dxux.*PI + dyuy.*PI +               2.*dxux.*(B + PI).*ut.^2).*ux.^5 - dxut.*(B + PI).*ut.*ux.^6 +            dxux.*(B + PI).*ux.^7)).*uy.^3 +          l.*(dxuy.*ut.^2.*(-B + PI.*(-3 + ut.^2)) -         dyut.*ut.*(B.*(8 - 8.*ut.^2 + ut.^4) + PI.*(6 - 5.*ut.^2 + ut.^4)).*ux +               (8.*B.*dxuy + 6.*B.*dyux + 6.*dxuy.*PI +            4.*dyux.*PI - (5.*B.*dxuy + 8.*B.*dyux + 4.*dxuy.*PI + 5.*dyux.*PI).*             ut.^2 + (dxuy + dyux).*(B + PI).*ut.^4).*ux.^2 +                       dyut.*ut.*(-5.*B - 3.*PI + 2.*(B + PI).*ut.^2).*          ux.^3 - (dxuy + dyux).*(-5.*B - 3.*PI + 2.*(B + PI).*ut.^2).*ux.^4 -         dyut.*(B + PI).*ut.*ux.^5 + (dxuy + dyux).*(B + PI).*ux.^6).*            uy.^4 +      l.*(dxut.*PI.*          ut.^3 + (8.*B.*dyuy +            6.*dyuy.*PI - (8.*B.*dyuy + (dxux + 5.*dyuy).*PI).*ut.^2 +            dyuy.*(B + PI).*ut.^4).*ux +                       dxut.*ut.*(-5.*B - 4.*PI + 2.*(B + PI).*ut.^2).*          ux.^2 + (5.*B.*dxux + 5.*B.*dyuy + 4.*dxux.*PI + 3.*dyuy.*PI -            2.*(dxux + dyuy).*(B + PI).*ut.^2).*ux.^3 -         2.*dxut.*(B + PI).*ut.*ux.^4 +               (2.*dxux + dyuy).*(B + PI).*ux.^5).*uy.^5 +      l.*((-dxuy).*PI.*ut.^2 +         dyut.*ut.*(-5.*B - 4.*PI + 2.*(B + PI).*ut.^2).*          ux - (dxuy + dyux).*(-5.*B - 4.*PI + 2.*(B + PI).*ut.^2).*ux.^2 -               2.*dyut.*(B + PI).*ut.*ux.^3 +         2.*(dxuy + dyux).*(B + PI).*ux.^4).*uy.^6 +      l.*ux.*(dyuy.*(5.*B + 4.*PI - 2.*(B + PI).*ut.^2) - dxut.*(B + PI).*ut.*ux +               (dxux + 2.*dyuy).*(B + PI).*ux.^2).*uy.^7 +      l.*(B + PI).*ux.*((-dyut).*ut + (dxuy + dyux).*ux).*uy.^8 +      dyuy.*l.*(B + PI).*ux.*uy.^9));
    
	
	
	Tlist(t) = max(max(T-0.7));
	uxlist(t) = max(max(ux));
	R = [R1 R2 R3 R4 R5];

	for tt=1:4
	
    [dt, fval, exitflag, output, jacobian] = fsolve(@eq, [DT Dux Duy DPI DB], optimset('Jacobian', 'on', 'DerivativeCheck', 'off', 'Display', 'off', 'Diagnostics', 'off'));

    DT  = dt(:,            1:  GRID_SIZE);
    Dux = dt(:,  GRID_SIZE+1:2.*GRID_SIZE);
    Duy = dt(:,2.*GRID_SIZE+1:3.*GRID_SIZE);
	DPI = dt(:,3.*GRID_SIZE+1:4.*GRID_SIZE);
	DB = dt(:,4.*GRID_SIZE+1:5.*GRID_SIZE);

    dtT = (DT-ux.*dxT-uy.*dyT)./ut;
	dtux = (Dux-ux.*dxux-uy.*dyux)./ut;
	dtuy = (Duy-ux.*dxuy-uy.*dyuy)./ut;
	dtPI = (DPI-ux.*dxPI-uy.*dyPI)./ut;
	dtB = (DB-ux.*dxB-uy.*dyB)./ut;
	
	
	
	T  = T0  + Dt * dtT;
	ux = ux0 + Dt * dtux;
    uy = uy0 + Dt * dtuy;
	ut = sqrt(1+ux.^2+uy.^2);
	PI = PI0 + Dt * dtPI;
	B = B0 + Dt * dtB;

    dxux = dx(ux);
    dyux = dy(ux);
    dxuy = dx(uy);
    dyuy = dy(uy);
	dxut = dx(ut);
	dyut = dy(ut);
    dxT = dx(T);
    dyT = dy(T);
	dxPI = dx(PI);
	dyPI = dy(PI);
	dxB = dx(B);
	dyB = dy(B);

    dxdxux = dxdx(ux);
    dxdyux = dx(dy(ux));
    dydyux = dydy(ux);
    dxdxuy = dxdx(uy);
    dxdyuy = dx(dy(uy));
    dydyuy = dydy(uy);
    dxdxT = dxdx(T);
    dxdyT = dx(dy(T));
    dydyT = dydy(T);



    phydro();
    pause(0.003);
  end
	end
















    phydro();
    %pause(0.033);

  function [M R] = eq(g)
    DT  = g(:,            1:  GRID_SIZE);
	Dux = g(:,  GRID_SIZE+1:2.*GRID_SIZE);
    Duy = g(:,  2.*GRID_SIZE+1:3.*GRID_SIZE);
    DPI = g(:,3.*GRID_SIZE+1:4.*GRID_SIZE);
	DB = g(:,4.*GRID_SIZE+1:5.*GRID_SIZE);

    

	L1 =(3.*DT.*T.^2.*(-1 + 3.*ut.^2) + 2.*Duy.*(B + PI).*ux + 6.*Dux.*T.^3.*ux +    2.*(Dux.*(B + PI) + 3.*Duy.*T.^3 + (DB + DPI).*ux).*uy)./(2.*ut);
	
    L2 = (1./(2.*ut.^4)).*(ut.^2.*(Duy.*(B + PI) + (B.*Duy + 2.*Duy.*PI + 3.*Dux.*T.^3).*        ux.^2 + 3.*T.^2.*ut.^2.*(Dux.*T + 3.*DT.*ux)) +       ((DB + DPI).*ut.^2 -       Dux.*(B + PI).*ux + (2.*B.*Dux + 4.*Dux.*PI + 3.*Duy.*T.^3).*ut.^2.*        ux + (DB + 2.*DPI).*ut.^2.*ux.^2 - Dux.*(B + 2.*PI).*ux.^3).*uy -       Duy.*(B + PI - 3.*PI.*ut.^2 + (B + 2.*PI).*ux.^2).*     uy.^2 + (DPI.*ut.^2 - Dux.*PI.*ux).*uy.^3 - Duy.*PI.*uy.^4);

    L3 = (1./(2.*ut.^4)).*(ut.^2.*(3.*Duy.*T.^3.*ut.^2 +       ux.*(DB + DPI + DB.*ux.^2)) + (9.*DT.*T.^2.*ut.^4 +       Duy.*PI.*(-1 + 2.*ut.^2).*ux).*uy +    ut.^2.*(3.*Duy.*T.^3 + (2.*DB + DPI).*ux).*uy.^2 -       Duy.*PI.*ux.*uy.^3 +    Dux.*(3.*T.^3.*ut.^2.*ux.*uy + PI.*(ut - ux).*(ut + ux).*(1 + uy.^2)) +       B.*((-Duy).*ux.*uy.*(1 - 4.*ut.^2 + ux.^2 + 2.*uy.^2) +       Dux.*((-ux.^2).*(1 + ux.^2 + 2.*uy.^2) + ut.^2.*(1 + 3.*ux.^2 + 2.*uy.^2))));
	
	L4 =(1./(4.*c.*T.^2.*ut.^2)).*(2.*B.*Duy.*h.*T.*ut.^2.*ux - 2.*B.*Duy.*l.*T.*ut.^2.*ux -    2.*Duy.*h.*PI.*T.*ut.^2.*ux - 2.*Duy.*l.*PI.*T.*ut.^2.*ux +    2.*B.*Duy.*h.*T.*ut.^4.*ux + B.*Duy.*l.*T.*ut.^4.*ux +       2.*Duy.*h.*PI.*T.*ut.^4.*ux + Duy.*l.*PI.*T.*ut.^4.*ux +    4.*B.*Duy.*h.*T.*ut.^2.*ux.^3 - 3.*B.*Duy.*l.*T.*ut.^2.*ux.^3 -    Duy.*l.*PI.*T.*ut.^2.*ux.^3 - 2.*B.*Duy.*h.*T.*ut.^4.*ux.^3 +    B.*Duy.*l.*T.*ut.^4.*ux.^3 -       2.*Duy.*h.*PI.*T.*ut.^4.*ux.^3 + 2.*B.*Duy.*h.*T.*ut.^2.*ux.^5 -    B.*Duy.*l.*T.*ut.^2.*ux.^5 + 2.*Duy.*h.*PI.*T.*ut.^2.*ux.^5 +       (2.*h.*ut.^2.*        ux.*((-(DB + 3.*DPI)).*T + (DB + DPI).*           T.*(ux.^4 - ut.^2.*(-1 + ux.^2)) +          2.*DT.*PI.*(3 - ux.^4 + ut.^2.*(-1 + ux.^2))) +                  Dux.*PI.*T.*(2.*           h.*(ut.^2.*(-3 + ut.^2) - (-1 + ut.^2).^2.*ux.^2 + (-4 + 3.*ut.^2).*              ux.^4 - 2.*ux.^6) +                          l.*(1 + ux.^2.*(5 + 8.*ux.^2 + 5.*ux.^4 + ux.^6 + ut.^4.*(3 + ux.^2) -                ut.^2.*(7 + 8.*ux.^2 + 2.*ux.^4)))) +                  B.*(4.*DT.*h.*ut.^2.*ux.*(-1 + ut.^2 - ux.^2).*(-1 + ux.^2) +          Dux.*T.*(l.*(1 +                ux.^2).*(1 + (3 - 3.*ut.^2 + ut.^4).*ux.^2 + (3 - 2.*ut.^2).*                 ux.^4 + ux.^6) -                                  2.*h.*(ut.^4.*(-1 + ux.^2) + ut.^2.*(1 + 2.*ux.^2 - 3.*ux.^4) +                2.*(ux.^4 + ux.^6))))).*uy +       Duy.*T.*     ux.*(PI.*(-2.*h + l -          2.*(h + l).*ut.^2 + (h.*(-8 + 6.*ut.^2) + l.*(6 - 5.*ut.^2 + ut.^4)).*           ux.^2 + (-4.*h + l.*(5 - 2.*ut.^2)).*ux.^4 + l.*ux.^6) +                  B.*(l.*(1 - ut.^4 + (-2 + ut.^2).^2.*ux.^2 - 2.*(-2 + ut.^2).*ux.^4 +             ux.^6) + 2.*h.*(3.*ut.^2.*(-1 + ux.^2) - 2.*(ux.^2 + ux.^4)))).*     uy.^2 +       (-2.*h.*(2.*DT.*(B + PI) - (DB + DPI).*T).*ut.^2.*ux.*(-1 + ux.^2) +       Dux.*T.*(2.*           h.*((-(B + PI)).*ut.^2 + (2.*B + (B + PI).*ut.^2).*ux.^2 -             2.*(B + PI).*ux.^4) +                          l.*(PI + (B + 3.*PI - 3.*PI.*ut.^2).*              ux.^2 + (3.*B + 5.*PI - 2.*(B + PI).*ut.^2).*ux.^4 +             2.*(B + PI).*ux.^6))).*uy.^3 +       Duy.*T.*     ux.*(4.*B.*h - B.*l + l.*PI +       2.*B.*l.*ut.^2 - (4.*B.*h - 3.*B.*l + 4.*h.*PI - 5.*l.*PI +          2.*l.*(B + PI).*ut.^2).*ux.^2 + 2.*l.*(B + PI).*ux.^4).*uy.^4 +       Dux.*l.*T.*ux.^2.*(-B + (B + PI).*ux.^2).*uy.^5 +    Duy.*l.*T.*ux.*(-B + (B + PI).*ux.^2).*uy.^6 +    2.*c.^2.*T.^4.*(Dux.*(ux - 2.*ut.^2.*ux + ux.^3) + Duy.*(1 + ux.^2).*uy));
	
	L5 = (1./(4.*c.*T.^2.*ut.^2)).*(T.*     ux.*(2.*c.^2.*Dux.*T.^3 + Duy.*PI.*(l + 2.*h.*ut.^2.*(-1 + ut.^2 - ux.^2)) +       B.*Duy.*(2.*h.*ut.^2.*(-3 + ut.^2 - ux.^2) + l.*(1 + ux.^2))) +       (2.*c.^2.*Duy.*T.^4.*(1 - 2.*ut.^2) +       Dux.*T.*ut.^2.*(-2.*B.*(h + l) +          2.*(h - l).*PI + (2.*h + l).*(B + PI).*ut.^2) +                  2.*h.*ut.^2.*(2.*DT.*(3.*B + PI) - (3.*DB + DPI).*           T + (-2.*DT.*(B + PI) + (DB + DPI).*T).*ut.^2).*ux -       Dux.*T.*(B.*(2.*h - l + 2.*(h + l).*ut.^2) +          PI.*(6.*h.*ut.^2 + l.*(-1 + ut.^4))).*              ux.^2 + 2.*h.*(2.*DT.*(B + PI) - (DB + DPI).*T).*ut.^2.*ux.^3 +       Dux.*T.*(4.*h.*PI + l.*(B + PI.*(-1 + 2.*ut.^2))).*ux.^4 -       Dux.*l.*PI.*T.*ux.^6).*uy +       T.*ux.*(2.*c.^2.*Dux.*T.^3 +       B.*Duy.*(-2.*h.*(-1 + ut.^2).^2 + 2.*h.*ut.^2.*ux.^2 +          l.*(5 - 7.*ut.^2 + 3.*ut.^4 - 3.*(-1 + ut.^2).*ux.^2)) +                  Duy.*PI.*(2.*h.*(2 + ut.^2).*(-ut.^2 + ux.^2) +          l.*(4 - 3.*ut.^2 + ut.^4 + ux.^2 - ux.^4))).*uy.^2 +       (2.*c.^2.*Duy.*T.^4 -       Dux.*T.*ut.^2.*(B.*l - 4.*h.*PI +          3.*l.*PI + (2.*B.*h + 2.*h.*PI - l.*PI).*ut.^2) +       2.*h.*(2.*DT.*(B + PI) - (DB + DPI).*T).*ut.^4.*ux +                  Dux.*T.*(-8.*B.*h + 6.*B.*l - 4.*h.*PI +          4.*l.*PI + (6.*B.*h - 5.*B.*l + 6.*h.*PI - 4.*l.*PI).*ut.^2 +          l.*(B + PI).*ut.^4).*ux.^2 +       2.*h.*(-2.*DT.*(B + PI) + (DB + DPI).*T).*ut.^2.*ux.^3 -                  Dux.*T.*(4.*B.*h - 5.*B.*l + 4.*h.*PI - 3.*l.*PI + 2.*l.*(B + PI).*ut.^2).*        ux.^4 + Dux.*l.*(B + PI).*T.*ux.^6).*uy.^3 +       Duy.*T.*     ux.*(B.*(h.*(6.*ut.^2 - 4.*(2 + ux.^2)) +          l.*(8 + ut.^4 + 5.*ux.^2 + ux.^4 - 2.*ut.^2.*(4 + ux.^2))) +                  PI.*(h.*(6.*ut.^2 - 4.*(1 + ux.^2)) +          l.*(6 + ut.^4 + 3.*ux.^2 + ux.^4 - ut.^2.*(5 + 2.*ux.^2)))).*uy.^4 +       (Dux.*(2.*B.*h + 2.*h.*PI - l.*PI).*T.*ut.^2 +       2.*h.*(-2.*DT.*(B + PI) + (DB + DPI).*T).*ut.^2.*ux -       Dux.*T.*(4.*B.*h - 5.*B.*l + 4.*h.*PI - 4.*l.*PI + 2.*l.*(B + PI).*ut.^2).*        ux.^2 +            2.*Dux.*l.*(B + PI).*T.*ux.^4).*uy.^5 +    Duy.*T.*ux.*(2.*PI.*(-2.*h + l.*(2 - ut.^2 + ux.^2)) +       B.*(-4.*h + l.*(5 - 2.*ut.^2 + 2.*ux.^2))).*uy.^6 +    Dux.*l.*(B + PI).*T.*ux.^2.*uy.^7 +       Duy.*l.*(B + PI).*T.*ux.*uy.^8);
	
	
	
    res = [L1 L2 L3 L4 L5] + R;

    % compute the gradients. hyhyhy



    DTL1 = (3.*T.^2.*(-1 + 3.*ut.^2))./(2.*ut);
    DTL2 =(9.*T.^2.*ux)./2;
	DTL3 =(9.*T.^2.*uy)./2;
	DTL4 =(h.*ux.*uy.*((-B).*(-1 + ux.^2).*(1 - ut.^2 + ux.^2 + uy.^2) +      PI.*(3 + ut.^2.*(-1 + ux.^2) + uy.^2 - ux.^2.*(ux.^2 + uy.^2))))./(c.*T.^2);
	DTL5 =(h.*ux.*uy.*(3.*B + PI - B.*ut.^2 - PI.*ut.^2 + B.*ux.^2 +      PI.*ux.^2 + (B + PI).*(ut - ux).*(ut + ux).*uy.^2 - (B + PI).*uy.^4))./(c.*    T.^2); 
    
	% ************************

    DuxL1 = (3.*T.^3.*ux + (B + PI).*uy)./ut;
    DuxL2 = (3.*T.^3.*ut.^2.*(ut.^2 + ux.^2) -    ux.*uy.*(B.*(1 - 2.*ut.^2 + ux.^2) +       PI.*(1 - 4.*ut.^2 + 2.*ux.^2 + uy.^2)))./(2.*ut.^4);
	DuxL3 = (3.*T.^3.*ut.^2.*ux.*uy + PI.*(ut - ux).*(ut + ux).*(1 + uy.^2) +    B.*((-ux.^2).*(1 + ux.^2 + 2.*uy.^2) + ut.^2.*(1 + 3.*ux.^2 + 2.*uy.^2)))./(2.*    ut.^4);
	DuxL4 = (1./(4.*c.*T.*ut.^2)).*(2.*c.^2.*T.^3.*ux.*(1 - 2.*ut.^2 + ux.^2) +       uy.*(B.*(-2.*           h.*(ut.^4.*(-1 + ux.^2) + ut.^2.*(1 + 2.*ux.^2 - 3.*ux.^4) +             2.*(ux.^4 + ux.^6)) + 2.*h.*(ut.^2 - 2.*ux.^2).*(-1 + ux.^2).*uy.^2 +                          l.*(1 + ux.^8 + 2.*ux.^6.*(2 - ut.^2 + uy.^2) +             ux.^2.*(4 - 3.*ut.^2 + ut.^4 + uy.^2 - uy.^4) +             ux.^4.*(6 + ut.^4 + 3.*uy.^2 + uy.^4 - ut.^2.*(5 + 2.*uy.^2)))) +                  PI.*(-2.*h.*(ux.^2 + ut.^4.*(-1 + ux.^2) + 2.*ux.^4.*(2 + ux.^2 + uy.^2) +             ut.^2.*(3 + uy.^2 - ux.^2.*(2 + 3.*ux.^2 + uy.^2))) +                          l.*(1 + uy.^2 +             ux.^2.*(5 + 8.*ux.^2 + ut.^4.*(3 + ux.^2) + 3.*uy.^2 +                ux.^2.*(ux.^2 + uy.^2).*(5 + ux.^2 + uy.^2) -                ut.^2.*(7 + 3.*uy.^2 + 2.*ux.^2.*(4 + ux.^2 + uy.^2)))))));
	DuxL5 = (1./(4.*c.*T.*      ut.^2)).*(ut.^2.*(-2.*B.*(h + l) +       2.*(h - l).*PI + (2.*h + l).*(B + PI).*ut.^2).*uy -       ux.^2.*(B.*(-l + 2.*(h + (h + l).*ut.^2) - l.*ux.^2) +       PI.*(6.*h.*ut.^2 - 4.*h.*ux.^2 + l.*(-1 + ux.^2 + (ut.^2 - ux.^2).^2))).*uy +       (ut.^2.*(4.*h.*PI -          l.*(B + 3.*PI) + (-2.*B.*h + (-2.*h + l).*PI).*           ut.^2) + (PI.*(l.*(-2 + ut.^2).^2 + h.*(-4 + 6.*ut.^2)) +          B.*(h.*(-8 + 6.*ut.^2) + l.*(6 - 5.*ut.^2 + ut.^4))).*ux.^2 -            (4.*B.*h - 5.*B.*l + 4.*h.*PI - 3.*l.*PI + 2.*l.*(B + PI).*ut.^2).*        ux.^4 + l.*(B + PI).*ux.^6).*uy.^3 +       ((2.*B.*h + 2.*h.*PI - l.*PI).*        ut.^2 - (4.*B.*h - 5.*B.*l + 4.*h.*PI - 4.*l.*PI + 2.*l.*(B + PI).*ut.^2).*        ux.^2 + 2.*l.*(B + PI).*ux.^4).*uy.^5 + l.*(B + PI).*ux.^2.*uy.^7 +       2.*c.^2.*T.^3.*ux.*(1 + uy.^2)); 	
    
	% ************************
    
    DuyL1 = (3.*T.^3.*ux + (B + PI).*uy)./ut;
    DuyL2 = (ut.^2.*(B + PI + (B + 2.*PI).*ux.^2) +    3.*T.^3.*ut.^2.*ux.*uy - (B + PI - 3.*PI.*ut.^2 + (B + 2.*PI).*ux.^2).*uy.^2 -    PI.*uy.^4)./(2.*ut.^4);
	DuyL3 = (3.*T.^3.*ut.^2.*(ut.^2 + uy.^2) -    ux.*uy.*(B + PI - 4.*B.*ut.^2 - 2.*PI.*ut.^2 +       B.*ux.^2 + (2.*B + PI).*uy.^2))./(2.*ut.^4);
	DuyL4 = (1./(4.*c.*T.*ut.^2)).*(ut.^2.*     ux.*(PI.*(l.*(-2 + ut.^2 - ux.^2) +          2.*h.*(-1 + ut.^2 - ut.^2.*ux.^2 + ux.^4)) +                  B.*(l.*(-2 + ut.^2 - ux.^2).*(1 + ux.^2) +          2.*h.*(1 + ut.^2 - (-2 + ut.^2).*ux.^2 + ux.^4))) +    2.*c.^2.*T.^3.*(1 + ux.^2).*uy +       ux.*(PI.*(-2.*h + l -          2.*(h + l).*ut.^2 + (h.*(-8 + 6.*ut.^2) + l.*(6 - 5.*ut.^2 + ut.^4)).*           ux.^2 + (-4.*h + l.*(5 - 2.*ut.^2)).*ux.^4 + l.*ux.^6) +                  B.*(l.*(1 - ut.^4 + (-2 + ut.^2).^2.*ux.^2 - 2.*(-2 + ut.^2).*ux.^4 +             ux.^6) + 2.*h.*(3.*ut.^2.*(-1 + ux.^2) - 2.*(ux.^2 + ux.^4)))).*     uy.^2 +       ux.*(4.*B.*h - B.*l + l.*PI +       2.*B.*l.*ut.^2 - (4.*B.*h - 3.*B.*l + 4.*h.*PI - 5.*l.*PI +          2.*l.*(B + PI).*ut.^2).*ux.^2 + 2.*l.*(B + PI).*ux.^4).*uy.^4 +    l.*ux.*(-B + (B + PI).*ux.^2).*uy.^6);
	DuyL5 = (1./(4.*c.*T.*ut.^2)).*(l.*PI.*     ux.*(1 + (4 - 3.*ut.^2 + ut.^4 + ux.^2 - ux.^4).*        uy.^2 + (6 + ut.^4 + 3.*ux.^2 + ux.^4 - ut.^2.*(5 + 2.*ux.^2)).*uy.^4 +       2.*(2 - ut.^2 + ux.^2).*uy.^6 + uy.^8) -       2.*((-c.^2).*T.^3.*uy.*(1 - 2.*ut.^2 + uy.^2) +       h.*PI.*ux.*(ut.^4.*(-1 + uy.^2) -          ut.^2.*(-1 + uy.^2).*(1 + ux.^2 + 3.*uy.^2) +          2.*(uy.^4 + uy.^6 + ux.^2.*uy.^2.*(-1 + uy.^2)))) +       B.*ux.*(-2.*        h.*(uy.^2 + ut.^4.*(-1 + uy.^2) + 2.*uy.^4.*(2 + ux.^2 + uy.^2) +          ut.^2.*(3 + ux.^2 - (2 + ux.^2).*uy.^2 - 3.*uy.^4)) +                  l.*(1 + (5 - 7.*ut.^2 + 3.*ut.^4).*uy.^2 + ux.^4.*uy.^4 +          uy.^4.*(8 + ut.^4 + 5.*uy.^2 + uy.^4 - 2.*ut.^2.*(4 + uy.^2)) +          ux.^2.*(1 + uy.^2.*(1 - ut.^2 + uy.^2).*(3 + 2.*uy.^2)))));
	
	
	% ************************
	
	DPIL1 = (ux.*uy)./ut;
	DPIL2 = (uy + 2.*ux.^2.*uy + uy.^3)./(2.*ut.^2);
	DPIL3 = (ux.*(1 + uy.^2))./(2.*ut.^2);
	DPIL4 = (h.*ux.*uy.*(-3 + ux.^4 - ut.^2.*(-1 + ux.^2) + (-1 + ux.^2).*uy.^2))./(2.*c.*T);
	DPIL5 = (h.*ux.*uy.*(-1 + uy.^2).*(1 - ut.^2 + ux.^2 + uy.^2))./(2.*c.*T);
	
	% ************************
	
	
	DBL1 =(ux.*uy)./ut;
	DBL2 =((1 + ux.^2).*uy)./(2.*ut.^2);
	DBL3 =(ux + ux.^3 + 2.*ux.*uy.^2)./(2.*ut.^2);
	DBL4 =(h.*ux.*(-1 + ux.^2).*uy.*(1 - ut.^2 + ux.^2 + uy.^2))./(2.*c.*T);
	DBL5 =(h.*ux.*uy.*(-3 + uy.^4 - ut.^2.*(-1 + uy.^2) + ux.^2.*(-1 + uy.^2)))./(2.*c.*T);

	% ********************
	
	M = sparse(zeros(5*GRID_SIZE*GRID_SIZE, 5*GRID_SIZE*GRID_SIZE));
	
	M = [diag((ut.*DTL1)(:)) diag((ut.*DuxL1)(:)) diag((ut.*DuyL1)(:)) diag((ut.*DPIL1)(:)) diag((ut.*DBL1)(:));
	diag((ut.*DTL2)(:)) diag((ut.*DuxL2)(:)) diag((ut.*DuyL2)(:)) diag((ut.*DPIL2)(:)) diag((ut.*DBL2)(:));
	diag((ut.*DTL3)(:)) diag((ut.*DuxL3)(:)) diag((ut.*DuyL3)(:)) diag((ut.*DPIL3)(:)) diag((ut.*DBL3)(:));
	diag((ut.*DTL4)(:)) diag((ut.*DuxL4)(:)) diag((ut.*DuyL4)(:)) diag((ut.*DPIL4)(:)) diag((ut.*DBL4)(:));
	diag((ut.*DTL5)(:)) diag((ut.*DuxL5)(:)) diag((ut.*DuyL5)(:)) diag((ut.*DPIL5)(:)) diag((ut.*DBL5)(:))];
	
	R=[(-R1-DTL1.*(ux.*dxT+uy.*dyT)-DuxL1.*(ux.*dxux+uy.*dyux)-DuyL1.*(ux.*dxuy+uy.*dyuy)-DPIL1.*(ux.*dxPI+uy.*dyPI)-DBL1.*(ux.*dxB+uy.*dyB))(:);
	(-R2-DTL2.*(ux.*dxT+uy.*dyT)-DuxL2.*(ux.*dxux+uy.*dyux)-DuyL2.*(ux.*dxuy+uy.*dyuy)-DPIL2.*(ux.*dxPI+uy.*dyPI)-DBL2.*(ux.*dxB+uy.*dyB))(:);
	(-R3-DTL3.*(ux.*dxT+uy.*dyT)-DuxL3.*(ux.*dxux+uy.*dyux)-DuyL3.*(ux.*dxuy+uy.*dyuy)-DPIL3.*(ux.*dxPI+uy.*dyPI)-DBL3.*(ux.*dxB+uy.*dyB))(:);
	(-R4-DTL4.*(ux.*dxT+uy.*dyT)-DuxL4.*(ux.*dxux+uy.*dyux)-DuyL4.*(ux.*dxuy+uy.*dyuy)-DPIL4.*(ux.*dxPI+uy.*dyPI)-DBL4.*(ux.*dxB+uy.*dyB))(:);
	(-R5-DTL5.*(ux.*dxT+uy.*dyT)-DuxL5.*(ux.*dxux+uy.*dyux)-DuyL5.*(ux.*dxuy+uy.*dyuy)-DPIL5.*(ux.*dxPI+uy.*dyPI)-DBL5.*(ux.*dxB+uy.*dyB))(:)]
	
	
	
	
    J = zeros(5*GRID_SIZE*GRID_SIZE, 5*GRID_SIZE*GRID_SIZE);
    for i=1:GRID_SIZE*GRID_SIZE
      x = rem(i, GRID_SIZE);
      y = mod(i, GRID_SIZE);

      du = sparse(zeros(GRID_SIZE, GRID_SIZE));
      du(i) = 1;
      
      Jxy = DTL1 .* du;
      J(                      1:  GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = DTL2 .* du;
      J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = DTL3 .* du;
	  J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DTL4 .* du;
	  J(3*GRID_SIZE*GRID_SIZE+1:4*GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DTL5 .* du;
	  J(4*GRID_SIZE*GRID_SIZE+1:5*GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  

      Jxy = DuxL1 .* du;
      J(                      1:  GRID_SIZE*GRID_SIZE,  GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = DuxL2 .* du;
	  J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,  GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DuxL3 .* du;
	  J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,  GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DuxL4 .* du;
	  J(3*GRID_SIZE*GRID_SIZE+1:4*GRID_SIZE*GRID_SIZE,  GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DuxL5 .* du;
	  J(4*GRID_SIZE*GRID_SIZE+1:5*GRID_SIZE*GRID_SIZE,  GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	 

      Jxy = DuyL1 .* du;
      J(                      1:  GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = DuyL2 .* du;
      J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = DuyL3 .* du;
      J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DuyL4 .* du;
	  J(3*GRID_SIZE*GRID_SIZE+1:4*GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DuyL5 .* du;
	  J(4*GRID_SIZE*GRID_SIZE+1:5*GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  
	  
	  Jxy = DPIL1 .* du;
      J(                      1:  GRID_SIZE*GRID_SIZE,3*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = DPIL2 .* du;
      J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,3*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = DPIL3 .* du;
      J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,3*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DPIL4 .* du;
	  J(3*GRID_SIZE*GRID_SIZE+1:4*GRID_SIZE*GRID_SIZE,3*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DPIL5 .* du;
	  J(4*GRID_SIZE*GRID_SIZE+1:5*GRID_SIZE*GRID_SIZE,3*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  
	  
	  Jxy = DBL1 .* du;
      J(                      1:  GRID_SIZE*GRID_SIZE,4*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = DBL2 .* du;
      J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,4*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = DBL3 .* du;
      J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,4*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DBL4 .* du;
	  J(3*GRID_SIZE*GRID_SIZE+1:4*GRID_SIZE*GRID_SIZE,4*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  Jxy = DBL5 .* du;
	  J(4*GRID_SIZE*GRID_SIZE+1:5*GRID_SIZE*GRID_SIZE,4*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
	  
    end
    % sqrt(sum(sum(res.*res)))
  end

  function r = phydro()
    frame_active = (mod(t-1,int16(1/Dt/SECONDS_PER_LIGHTYEAR/30)) == 0);
    if VIDEO & ~frame_active 
      return
    end
    if VIDEO
      disp(t);
    end
    plotsx = 4;
    plotsy = 3;

    subplot(plotsy, plotsx, 1);
    image((T-TEMP_OFFSET) * 256  + 128); 
    
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vt,frame);
    % end

    subplot(plotsy, plotsx, 2);
    image(ux*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvx,frame);
    % end

    subplot(plotsy, plotsx, 3);
    image(uy*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvy,frame);
    % end

    subplot(plotsy, plotsx, 4);
    image(sqrt(1+ux.*ux + uy.*uy)*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vv,frame);
    % end

    subplot(plotsy, plotsx, 5);
    image(dtT*4+128);

    subplot(plotsy, plotsx, 6);
    image(dtux*4+128);
    subplot(plotsy, plotsx, 7);
    image(dtuy*4+128);

    %truesize(fig);
    if (VIDEO) & (frame_active)
      frame = getframe(gcf);
      writeVideo(vt,frame);
    end

    
    % subplot(plotsy, plotsx, 9);
    % image(dtdtvx*8 + 128);
    % subplot(plotsy, plotsx, 10);
    % image(dtdtvy*8 + 128);
    % subplot(plotsy, plotsx, 11);
    % image(Pi*2 + 128);
    % subplot(plotsy, plotsx, 12);
    % image(B*2 + 128);
    return;

  end

  % [-1 0 1
  function r = dx(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = - arg(GRID_SIZE,y)/2 + arg(2, y)/2;
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  -arg(GRID_SIZE-1,y)/2 + arg(1, y)/2;
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = (arg(x+1,y) - arg(x-1,y))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dy(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = - arg(x,GRID_SIZE)/2 + arg(x,2)/2;
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  -arg(x,GRID_SIZE-1)/2 + arg(x,1)/2;
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = (arg(x,y+1) - arg(x,y-1))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dx2(arg)
    r = d1 * arg;
  end
  function r = dy2(arg)
    r = arg * (d1.');
  end
  function r = dxdx2(arg)
    r = d2 * arg;
  end
  function r = dydy2(arg)
    r = arg * (d2.');
  end
  function r = dxdx(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = arg(GRID_SIZE,y) + arg(2, y) - 2*arg(1,y);
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  arg(GRID_SIZE-1,y) + arg(1, y) - 2*arg(GRID_SIZE,y);
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = arg(x+1,y) + arg(x-1,y) - 2*arg(x,y);
      end
    end 
    r = r / DX / DX;
  end
  function r = dydy(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = arg(x,GRID_SIZE) + arg(x,2) - 2*arg(x,1);
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  arg(x,GRID_SIZE-1) + arg(x,1) - 2*arg(x,GRID_SIZE);
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = arg(x,y+1) + arg(x,y-1) - 2*arg(x,y);   
      end
    end 
    r = r / DX / DX;
  end

end


